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Abstract — We  report  a novel  method  for  accurately  com- 
puting the  modes  of  an  arbitrarily-shaped  hollow  wave- 
guide. Our  method  uses  a point-based  (Nystrom)  dis- 
cretization of  an  integral  operator  over  the  waveguide 
aperture  to  compute  the  modes. 

I.  Introduction 

This  paper  describes  a technique  for  numerically  com- 
puting the  modes  of  an  arbitrarily-shaped  cylindri- 
cal waveguide  with  perfectly-conducting  walls  and  an 
isotropic,  homogeneous  core.  The  method  is  different 
from  standard  methods  in  that  we  formulate  the  prob- 
lem in  terms  of  the  eigenfunctions  of  an  integral  operator 
over  the  waveguide  aperture  and  it  is  unique  in  that  we 
employ  a high-order,  point-based  (Nystrom)  discretiza- 
tion to  obtain  numerical  solutions. 

The  Nystrom  method  is  a method  for  solving  inte- 
gral equations.  In  contrast  to  a method  of  moments  dis- 
cretization, a Nystrom  discretization  of  a function  on  a 
surface  S is  simply  a tabulation  of  function  values  at  a 
discrete  set  of  points  on  S . Integrals  are  approximated 
by  weighted  sums  of  function  evaluations.  Specifically, 
we  approximate  the  integral  of  a function  / ( x ) as 

N 

dx  f (x)  ~ £“>./(*.),  (i) 

i-i 

where  Xi  is  the  ith  abscissa  of  an  iV-point,  high-order 
quadrature  rule  and  u>i  is  the  associated  quadrature 
weight. 

The  conventional  Nystrom  method  is  a simple  and  ef- 
ficient method  for  solving  integral  equations  with  non- 
singular kernels.  When  the  integral  kernel  is  singular 
(as  is  generally  the  case  for  Green  functions),  one  needs 
to  introduce  local  corrections  in  order  to  compensate  for 
the  fact  that  a quadrature  rule  for  regular  functions  can- 
not integrate  singular  functions  with  high-order  accuracy. 
Further  details  regarding  local  corrections  and  their  use 
in  the  Nystrom  method  for  solving  scattering  problems 
may  be  found  in  [1]. 


The  modes  of  an  isotropic,  homogeneous-core  cylin- 
drical waveguide  with  perfectly  conducting  walls  can  be 
classified  [2]  into  three  categories  — transverse  magnetic 
(TM),  transverse  electric  (TE),  and  transverse  electro- 
magnetic (TEM).  Modes  in  each  category  are  related  to 
the  modes  of  a scalar  waveguide  problem.  The  TM  (TE) 
modes  are  derivable  from  the  scalar  modes  of  the  same 
waveguide  assuming  Dirichlet  (Neumann)  boundary  con- 
ditions on  its  walls.  The  TEM  modes  are  derivable  from 
solutions  to  the  2d  Laplace  equation  on  a cross  section  of 
the  guide. 

The  arrangement  of  the  remainder  of  the  paper  is  as 
follows:  section  II  describes  how  to  obtain  eigenmodes 
and  eigenvalues  for  the  scalar  waveguide  problem  with 
Dirichlet  or  Neumann  boundary  conditions;  section  III 
describes  how  to  solve  the  2d  Laplace  equation  inside  the 
waveguide  aperture;  section  IV  describes  a method  for 
computing  the  vector  waveguide  modes  from  the  scalar 
traveling  modes  and  2d  electrostatic  modes;  and,  finally, 
in  section  V we  present  some  results  obtained  from  a soft- 
ware implementation  of  these  methods. 

II.  Scalar  waveguide  modes 

The  scalar  waveguide  modes  un  (x)  satisfy  the  2d  scalar 
wave  equation 

(Vi  + A;2  - Pi)  un  (x)  = 0 (2) 

inside  the  waveguide  aperture  W and  the  correct  bound- 
ary conditions  on  its  boundary  d W.  They  are  also  eigen- 
functions of  the  H operator  [3]  defined  as 

H (x,  x')  = £ ~[£Un  (x) Un  (x0  • (3) 

n 

In  these  equations, 


is  the  modal  impedance,  k is  the  free-space  propagation 
constant,  and  f3n  is  the  propagation  constant  for  the  nth 
waveguide  mode. 
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Using  the  fact  that  the  modes  form  a complete  and 
orthonormal  set  of  real  functions  on  W , one  can  show 
that  the  function  G defined  by 

G(x,  x")  = f ds'  H H (x' ,x") , (5) 

Jw 

is  the  Green  function  appropriate  to  the  inside  of 
the  waveguide  because  it  obeys  the  wave  equation 
(Vj_  + A;2)  G (x,x")  = — <5(2)  (x,x")  inside  W and  the 
correct  boundary  conditions  on  dW.  Furthermore,  the 
eigenfunctions  of  the  integral  operator 

J ds " G (x,  x")  = (7^)  Un  W Jw  ds>  un  (x') 

(6) 

are  the  modes  of  the  waveguide  and  the  eigenvalue  cor- 
responding to  the  nth  mode  un  (x)  is 


Therefore,  our  procedure  for  computing  the  scalar  trav- 
eling modes  of  a waveguide  will  consist  of  first  comput- 
ing a discretized  representation  of  the  integral  opera- 
tor  fw  dsn  G (x,  x")  and  then  using  a numerical  eigen- 
value routine  to  determine  discretized  representations  (of 
a finite  set)  of  the  modes  and  the  corresponding  prop- 
agation constants.  The  remainder  of  this  sections  de- 
scribes a method  for  obtaining  discretized  representations 
of  fw  ds"  G (x,  x")  with  Dirichlet  or  Neumann  boundary 
conditions  on  dW . 

Start  with  a solution  to  the  inhomogeneous  (trans- 
verse) wave  equation  -fi  /c2)  G (x,x')  = — 6^  (x,  x'). 
We  will  use 

G(x,x')  = -~Y0{k  |x-x'|),  (8) 

where  To  is  the  second  kind  Bessel  function  of  order  zero. 
To  this  solution  we  can  always  add  solutions  F (x,x')  to 
the  homogeneous  wave  equation  (V2  + k 2)  F (x,  x')  = 0. 
Our  objective  is  to  find  a solution  F that  makes 

G (x,  x')  = G (x,  x')  + F (x,  x')  (9) 

obey  the  boundary  conditions  on  dW.  Unlike  the  Green 
function  for  2d  scattering  in  an  unbounded  region,  this 
Green  function  is  real  valued. 

A.  Dirichlet  Case 

Since  we  care  only  about  G inside  W,  we  can  arrange 
any  distribution  of  charges  cr  outside  of  W to  make  G 
obey  the  boundary  conditions  on  dW . The  simplest  so- 
lution is  to  put  them  on  an  artificial  boundary  T that  is 
outside  dW  by  an  infinitesimal  distance.  Then 

F(x,x')  = J dl"  G(x,x")cr(x",x/)  (10) 


obeys 

(V^  + k2)  F (x,  x') 

= J dl"  (Vl  + k 2)  G (x,  x")  a (x",  x')  = 0 (11) 

since  (V2  + k 2)  G (x,  x")  = 0 for  all  x € W and  x;/  G T. 
Taking  the  limit  as  T — * dW  we  get 

G (x,  x;)  = G (x,  x;)  + f dl"  G(x,x")<j(x",x!)  . 

Jew 

(12) 


A simple  interpretation  of  this  equation  is  as  follows: 
cr  (x"  e dW)  x;  e W)  is  the  charge  distribution  induced 
on  the  walls  of  the  waveguide  by  a unit  charge  at  x'  G W\ 
the  total  potential  at  x G dW  is  the  sum  of  the  potential 
from  the  original  unit  charge,  namely  G (x,x'),  and  the 
potential  produced  by  the  induced  charge  distribution  on 
the  waveguide  walls,  namely 


/ . 


dl"  G(x,x'V(x",x'). 


(13) 


In  the  Dirichlet  case,  the  total  potential  must  vanish 
everywhere  on  dW . To  enforce  this  condition,  we  will 
demand  that  the  inner  product  of  the  potential  on  dW 
with  each  function  fk  (x  G dW)  from  a suitable  set  of 
testing  functions  must  vanish,  i.e., 


0 = / dl  fk  (x)G(Xjx') 

Jaw 

+ / dl  fk  (x)  [ dl"  G(x,x")cr(x//,x/) 
JdW  JdW 

for  every  point  x'  G W. 

We  can  write  this  condition  in  matrix  form  as 


(14) 


o = fdWfldW  (^QdW,W  QdW.dW^dW^dW.W^  ^5) 


where  %dw'w  represents  the  discretized  form  of  cr  and 
GdWtdW  and  GdW,W 

are  discretized  representations  of 
the  kernel  G with  local  corrections  [1].  Ddw  is  a diagonal 
matrix  of  quadrature  weights  for  integrals  over  dW.  [In 
general,  our  notation  involving  W and  dW  superscripts 
is  meant  to  indicate  the  domain  (s)  of  the  coordinate  vari- 
able (s).  For  diagonal  quadrature  weight  matrices  fi,  only 
one  superscript  is  used  for  notational  compactness  with 
the  understanding  that  the  two  domains  are  always  the 
same.]  Since  fdW  is  arbitrary,  the  solution  to  (15)  is 


^dw,w  _ _ ^Qdw,dw GdW}W  (16) 


which  means  that  the  discretized  form  of  fw  ds"  G (x,  x") 
is 


Gw'wVtw  — Gw,dw  (pdwjwy1  qQW.WqW  ^7) 


86 


ACES  JOURNAL,  VOL.  17,  NO.  1,  MARCH  2002,  SI:  APPROACHES  TO  BETTER  ACCURACY/RESOLUTION  IN  CEM 


B.  Neumann  Case 

As  in  the  Dirichlet  case,  we  arrange  sources  on  T in 
just  the  right  amount  to  make  the  total  potential  obey 
the  boundary  conditions  on  dW . In  the  Neumann  case, 
however,  a dipole  distribution  ji  has  some  advantages  over 
a charge  distribution1. 

The  potential  at  any  point  x EW  due  to  a dipole  dis- 
tribution /i  on  F is 

E(x,x')  = J dl"  (e"  • V"G(x,x")) /z(x",x') , (18) 

where  e"  is  the  unit  normal  to  dW  at  x"  in  the  plane  of, 
but  pointing  away  from,  W.  For  the  same  reason  as  be- 
fore,  F obeys  the  homogeneous  wave  equation  whenever 
x EW.  The  edge  normal  derivative  of  F at  x gF  is 

e-  VF(x,x')  = (e  • V)  J dl"  (e"  • V"G(x,x"))m(x",x') 

(19) 

The  normal  derivative  of  a dipole  layer  potential  is  con- 
tinuous across  the  boundary  so  we  can  take  the  limit 
F — > d IT,  replacing  T by  dW  in  the  above  expression. 

In  the  Neumann  case,  the  edge  normal  derivative  of 
the  total  potential  must  vanish  everywhere  on  dW.  We 
will  demand  that  the  inner  product  of  the  edge  normal 
derivative  of  the  potential  on  dW  with  each  function  from 
a suitable  set  of  testing  functions  must  vanish.  In  other 
words,  for  each  testing  function  fk  (x  G dW),  we  require 
that 


derivative  operators  with  local  corrections.  Qw  and  £ldW 
are  diagonal  matrices  of  quadrature  weights  for  integrals 
over  W and  dW,  respectively.  The  solution  to  (21)  is 

Mdw'w 

= - ([(e  • V)  (e"  • V'G)]dW'dW  n9lv)_1  (e  • VG)dw’w , 

(22) 

which  means  that  the  discretized  form  of  Jw  ds"  G (x,  x") 
in  the  Neumann  case  is 

GwWsiw-{e-  VG)dW'w 

([(e  • V)  (e"  • V'G)]dW'dwyl  (e  ■ VG)dW'w  Q.w . (23) 

III.  Solving  Laplace’s  equation  in  W 

When  the  waveguide  aperture  is  multiply-connected, 
there  exist  non-trivial  solutions  to  the  2d  Laplace  equa- 
tion in  W.  Such  solutions  correspond  to  the  TEM  modes 
in  the  waveguide. 

What  follows  is  a procedure  for  computing  the  solution 
to  the  2d  Laplace  equation  in  W for  an  arbitrary  potential 
distribution  on  the  boundaries.  The  same  procedure  has 
an  obvious  extension  to  3 d which  could  be  used  to  solve 
electrostatic  problems  inside  multiply-connected  cavities. 

For  a given  boundary  value  function  b(xEdW))  we 
desire  to  find  a function  0 (x  EW)  that  satisfies 

0 (x)  = b (x)  for  x edW)  and  (24) 

V20  (x)  =0  for  x eW.  (25) 


0 = [ dl  fk  (x)  (e  • VG  (x,  x'))  + 

JdW 


The  solution  for  0 (x)  can  be  written  as  a 2d  single-layer 
potential  with  an  unknown  source  distribution 


dlfk(x)  (8-V) 


V"G(x,x"))m(x",x') 

(20) 


for  every  point  x'  G W. 

We  can  write  this  in  matrix  form  as 

0 = f9Wn9W  ((e  • VG)dw'w  + 


0(x)=  f dl!  log  |x  — x;|  a (x') . (26) 

JdW 

This  potential  automatically  satisfies  the  second  condi- 
tion above  because  V2log|x  — x'|  = 0 for  x EW  and 
xfEdW.  The  single-layer  potential  density  a (x')  is  deter- 
mined by  the  condition  that  0 (x)  = b (x)  on  the  bound- 
ary, i.e., 


[(e  • V)  (e"  ■ v"G)]9W'dW  SldwMdw'w^  (21) 


I ■ 

JdW 


dl'  log  |x  — x'|  u (x;)  = 6 (x) . 


(27) 


where  Mdw'w  represents  the  discretized  form  of  /z, 
and  [(©  • V)  (e"  ■ v"G))dW'dW  and  (e  • VG)dw'w  are  dis- 
cretized representations  of  the  corresponding  edge  normal 

^he  problem  with  using  a charge  distribution  in  the  Neumann 
case  is  that  it  diverges  at  acute  angle  corners,  such  as  in  a square 
waveguide.  A Id  quadrature  rule  designed  to  integrate  regular  func- 
tions will  not  be  high  order  for  such  a charge  distribution.  For  wave- 
guide apertures  with  smooth  edges,  such  as  a circular  guide,  this 
is  not  a problem  and  using  a charge  distribution  may  well  better 
since  the  integral  operator  is  second  kind. 

When  using  a charge  distribution  cr  for  the  Neumann  case,  beware 
of  the  fact  that  the  limit  of  fr  dl"  (e  • VG  (x,  x"))  cr  (x",  x')  as 
P — * dW  is  (x,  x;)  + fQW  dl"  (e  • VG  (x,  x"))  a (x",  x') . 


If  ipw^w  and  Lw'dw  are  the  discretized  representa- 
tions of  logjx  — x'l  (with  local  corrections)  for  x G dW 
and  x G W,  respectively,  and  EdW , and  BdW  are 

the  discretized  representations  of  0(x),  a (x)  and  6(x), 
respectively,  then 

j^QW.dW Q^dW  Y^dW  ___  gdW  ^8) 

and  the  solution  for  becomes 

tyW  __  jJV.dW q3W YpW  _ ^W,dW  ^dW,dW^~l  gdW 

(29) 
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To  get  TEM  modes,  W must  be  multiply  connected, 
i.e., 


dW  ^dWiUdW2...U  dWn  (30) 


solution  to  the  linear  system 


£ A (*<,*,)  Ate)  = d*k-^ . 


(38) 


where  the  dWi  are  unconnected  boundaries  and  n > 2. 
There  is  a TEM  mode  corresponding  to  each  of  the  n — 1 
independent  boundary  functions  described  by 


bik  (x)  = 0 for  all  k except  k = i (31) 

bu  (x)  = 1 (32) 

for  i = 1, 2, . . . ,n  — 1. 


IV.  Vector  waveguide  modes 

The  transverse  electric  components  of  the  TM,  TE, 
and  TEM  vector  waveguide  modes  are  derivable  from  the 
scalar  waveguide  modes  according  to 


u 


TM 


U 


TE 


(X) 


(x) 


V</>„ 


n X Vlpn 


nlEM  (x)  OC  VCn, 


(33) 

(34) 

(35) 


where  ipn)  ipni  and  ( n are  the  scalar  Dirichlet,  Neumann, 
and  Laplace  modes,  respectively.  The  transverse  mag- 
netic component  [2]  is 


Hi  = ±Z~1n  x Ej_,  (36) 


where 


for  n E TM  modes 
for  n E TEM  modes  (37) 
for  n E TE  modes 

is  the  modal  impedance  and  e and  \i  are  the  dielectric 
constant  and  magnetic  permeability,  respectively. 

We  need  discretized  representations  of  the  surface  gra- 
dient to  effect  the  transformation  from  scalar  modes  to 
vector  modes  given  above.  This  section  shows  how  to 
represent  the  surface  gradient  operator  (on  an  arbitrary 
surface)  in  matrix  form.  Left  multiplying  a matrix  repre- 
senting a scalar  surface  function  by  the  matrix  represent- 
ing the  surface  gradient  produces  a discretized  represen- 
tation of  the  surface  gradient  of  the  scalar  function. 

First  consider  the  linear  derivative  operator.  In  the 
spirit  of  the  high-order  Nystrom  method,  we  will  de- 
mand that  the  discretized  derivative  operator  return  ex- 
act results  at  a particular  set  of  sample  points  for  each 
function  in  a set  of  suitable  functions.  In  other  words, 
if  we  are  given  a set  of  points  Xj  on  a curve  (7,  then 
A ij  = A (xi,xj)  is  a high-order  discretized  representa- 
tion of  the  differential  operator  ~ at  Xi  on  C if  it  is  the 


for  suitable  testing  functions  fk  (a:). 

We  can  also  make  a connection  with  the  Nystrom 
method  by  re-expressing  the  derivative  as  an  integral  op- 
erator and  applying  the  standard  procedure  [1]  for  com- 
puting local  corrections.  If  we  write  the  linear  derivative 
of  / (z)  on  the  curve  C as 

^ = ljcdx'6{x-x')f{x,)’  w 

then  the  discretized  representation  of  the  differential  op- 
erator ^ on  C is  obtained  by  solving  the  linear  system 


A (xi,Xj)  fk  (xj) 

3 


L 


dx ' S(x  — x;)  fk  (xf) 


dfk  (gi) 
dx 


(40) 


for  A (xiyXj)  using  suitable  testing  functions  fk  (x). 
Clearly,  A and  A are  related  by  Aft  = A.  The  only 
difference  between  this  linear  system  and  the  local  cor- 
rection linear  system  in  [1]  is  that  computing  the  right 
hand  side  of  (40)  only  requires  evaluating  derivatives  of 
the  testing  functions  at  the  sample  points  instead  of  eval- 
uating inner  products  of  the  kernel  with  testing  functions. 

We  have  encountered  operators  similar  to  this  in  scat- 
tering problems  before.  The  hypersingular  operators 
(n  • V)  f ds*  (n'-V')  G (x^xt)  and  (n  x V)  f dsf  (n'xV') 
G (x,  a;'),  which  appear  in  boundary  integral  formulations 
of  scalar  and  electromagnetic  scattering,  respectively,  are 
pseudo-differential  operators.  Like  these  operators,  the 
discretized  gradient  operator  must  be  used  with  extreme 
caution  (and  avoided  whenever  possible)  because  it  tends 
to  amplify  rather  than  attenuate  numerical  “noise” . 

Now  consider  surface  derivative  operators.  Let  ti  (x) 
and  t2  (x)  be  independent  unit  tangent  vectors  on  the 
surface  S.  By  analogy  with  the  linear  derivative  opera- 
tor, we  obtain  a locally  corrected  matrix  representation 
Am  (xi,Xj)  of  the  surface  gradient  operator  ■ V on  S 
by  solving  the  linear  system 


^ ^ A/x  (xj, Xj)  fk  (xj ) — (xj)  • V fk  (xj)  (41) 
j 


with  p = 1,2  using  suitable  testing  functions  fk  (x).  If 
these  testing  functions  afford  a high-order  approximation 
to  a scalar  function  ip  (x)  on  5,  then  the  vector  A^ 
represents  a high-order  approximation  to  (x)  • V'tp  (x), 
with  (AM){j  = Am  (xi,xj)  and  = ip  (x<) . 

The  matrix  representation  of  n x V is  obtained  by  re- 
placing tM  (xj)  • V with  (xj)  • (n  x V)  in  (41). 
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TABLE  I TABLE  II 

Computed  and  Exact  Eigenvalues  (/ 3m,n/k ) for  a 1.1  A x Computed  and  Exact  Eigenvalues  (/3m,n/fc)  for  a 1.1A  x 

0.75A  Waveguide  Satisfying  Dirichlet  Boundary  Conditions  0.75A  Waveguide  Satisfying  Neumann  Boundary  Conditions 


4x4 

6x6 

8x8 

Exact 

m 

n 

0.592201 

0.590741 

0.590716 

0.590715 

1 

1 

0.512091i 

0.520366* 

0.520471* 

0.520472* 

2 

1 

0.969128* 

0.991772* 

0.992159* 

0.992164* 

1 

2 

1.219852* 

1.159101* 

1.142176* 

1.141906* 

3 

1 

1.629548i 

1.265721* 

1.266571* 

1.266580* 

2 

2 

1.997986* 

1.574208* 

1.624191* 

1.623971* 

3 

2 

2.080501i 

1.637739* 

1.648007* 

1.658382* 

4 

1 

2.309772* 

1.817060* 

1.791106* 

1.790701* 

1 

4 

2.546563* 

1.929628* 

1.956519* 

1.956130* 

2 

4 

2.694088* 

1.981966* 

2.010146* 

2.020783* 

4 

2 

3.010537* 

2.237905* 

2.204944* 

2.147029* 

5 

1 

3.114740* 

2.385630* 

2.279746* 

2.204428* 

3 

3 

3.321089* 

2.444649* 

2.497769* 

2.437841* 

5 

2 

3.469637* 

2.503110* 

2.501035* 

2.511132* 

4 

3 

3.642493* 

2.696163* 

2.537604* 

2.513508* 

1 

4 

3.744043* 

2.836216* 

2.560426* 

2.633924* 

2 

4 

V.  Results 

We  implemented  these  techniques  in  software  to  com- 
pute the  modes  of  arbitrarily-shaped,  closed  waveguides. 
Our  code  requires  two  inputs.  The  first  is  a description 
of  the  waveguide  aperture.  The  aperture  is  described  in 
terms  of  one  or  more  quadrilateral  or  triangular  patches. 
The  mesh  always  covered  the  planar  waveguide  aperture 
exactly  in  order  to  preserve  the  ability  to  achieve  high- 
order  convergence  in  the  solution.  We  locate  discretiza- 
tion points  on  each  of  these  patches  according  to  a high- 
order  2d  quadrature  rule.  Discretization  points  on  the 
boundary  are  located  according  to  a Id  quadrature  rule 
of  the  same  order.  The  number  of  such  points  is  deter- 
mined by  the  second  input,  the  order  of  the  quadrature 
rule.  This  value  also  determines  the  maximum  order  of 
the  testing  functions  used  to  compute  local  corrections. 
The  output  of  the  code  consists  of  numerically  computed 
eigenmodes  and  eigenvalues. 

We  have  tested  the  code  by  using  it  to  compute  modes 
of  several  waveguides  with  simple  cross  sections.  Sample 
results  for  two  waveguides  are  presented  in  this  section. 
The  first  is  a rectangular  waveguide,  a problem  for  which 
analytical  solutions  are  available.  We  compare  computed 
propagation  constants  (eigenvalues)  to  analytical  results. 
The  second  problem  is  a rectangular  waveguide  contain- 
ing two  square  conductors.  Such  a waveguide  has  two 
TEM  modes  in  addition  to  its  TM  and  TE  modes.  We 
list  computed  propagation  constants  and  plot  the  low- 
est modes.  Similar  results  have  been  obtained  on  other 
waveguide  shapes  including  circular  waveguide,  circular 
coaxial  waveguide,  and  rectangular  waveguide  with  sep- 
tum. 


4x4 

6x6 

8x8 

Exact 

m 

n 

0.889957 

0.890835 

0.890726 

0.890724 

1 

0 

0.745280 

0.745429 

0.745358 

0.745356 

0 

1 

0.587578 

0.590805 

0.590715 

0.590715 

1 

1 

0.443112 

0.412701 

0.416826 

0.416598 

2 

0 

0.507053* 

0.527063* 

0.520106* 

0.520472* 

2 

1 

0.758857* 

0.883037* 

0.881869* 

0.881917* 

0 

2 

0.919138* 

0.922549* 

0.925573* 

0.927094* 

3 

0 

0.928154* 

0.994163* 

0.992079* 

0.992164* 

1 

2 

1.028942* 

1.141903* 

1.139574* 

1.141906* 

3 

1 

1.288202* 

1.271016* 

1.266358* 

1.266580* 

2 

2 

1.470524* 

1.623864* 

1.536828* 

1.518482* 

4 

0 

1.574991* 

1.713198* 

1.622186* 

1.623971* 

3 

2 

1.606358* 

1.739490* 

1.681116* 

1.658382* 

4 

1 

1.769745* 

1.772999* 

1.731271* 

1.732051* 

0 

3 

1.894675* 

1.857114* 

1.789719* 

1.790701* 

1 

3 

1.932489* 

1.955014* 

1.956130* 

2 

3 

Tables  I and  II  list  eigenvalues  for  an  a x b rectangu- 
lar waveguide  with  a = 1.1  A and  b = 0.75 A.  For  com- 
putational purposes  the  waveguide  aperture  was  defined 
by  a single  rectangular  patch.  The  first  column  of  each 
table  gives  the  eigenvalues  computed  using  a 16-point, 
high-order  quadrature  rule,  namely  a product  rule  con- 
structed from  two  Id  Gauss-Legendre  rules  each  using  4 
points.  The  second  and  third  columns  give  the  16  low- 
est eigenvalues  computed  using  high-order  product  rules 
involving  36  and  64  points,  respectively. 

The  fourth  column  gives  the  exact  values  of  the  16  low- 
est eigenvalues.  The  analytical  solutions  [2]  for  a rectan- 
gular waveguide  are  well  known.  The  modes  satisfying 
Dirichlet  boundary  conditions  take  the  form 

<Pmn  (%,  y ) = sin  (m7rf ) sin  (r*7r|)  for  to, n = 1, 2, 3, . . . 

(42) 

Similarly,  the  Neumann  modes  are 
'ij?mn  (^>  y)  = cos  (m7r-)  cos  (titt |)  for  m,  n = 0, 1, 2, . . . 

(43) 

Accordingly,  the  exact  eigenvalues  for  both  Dirichlet  and 
Neumann  boundary  conditions  are 

¥=\A <44> 

The  last  two  columns  show  the  values  of  m and  n for  each 
computed  mode.  Four  Neumann  modes  and  one  Dirichlet 
mode  are  propagating  modes,  the  rest  are  evanescent. 

Not  surprisingly,  the  lower  the  mode,  the  more  accu- 
rate the  eigenvalue  for  a given  quadrature  rule.  With  a 64 
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TABLE  III 

Computed  Eigenvalues  (P/k)  vs.  Discretization  for  a 5A  x 
3 A Waveguide  Containing  Two  Interior  Conductors 


First  10  TM  modes 
4x4  6x6 

First  10  TE  modes 
4x4  6x6 

0.899770 

0.898954 

0.996147 

0.996122 

0.889666 

0.888867 

0.991385 

0.991335 

0.880594 

0.879772 

0.986280 

0.986221 

0.877820 

0.877004 

0.985071 

0.985030 

0.866849 

0.866051 

0.968195 

0.968162 

0.858987 

0.858275 

0.957283 

0.957146 

0.824123 

0.822946 

0.948050 

0.947774 

0.814496 

0.813010 

0.946083 

0.946051 

point  quadrature  rule,  the  lowest  eigenvalues  are  accurate 
to  almost  6 digits.  The  higher  the  spatial  frequency  con- 
tent of  the  mode,  however,  the  lower  the  accuracy  of  the 
computed  result.  This  is  evident  in  the  list  of  Dirichlet 
eigenvalues  computed  using  a 16-point  quadrature  rule. 
It  also  holds  in  both  cases  for  the  higher  modes  that  were 
computed2  but  are  not  shown  in  the  tables. 

The  second  sample  problem  is  a 3A  x 5A  rectangular 
waveguide  containing  two  1A  x 1A  square  conductors.  Ta- 
ble III  lists  the  computed  eigenvalues  for  the  first  ten 
Dirichlet/TM  and  Neumann/TE  modes  of  the  guide.  All 
computations  were  performed  using  a mesh  consisting  of 
thirteen  1A  x 1A  patches  arranged  on  a Cartesian  grid. 
The  columns  labeled  ‘4x4’  and  ‘6x6’  give  results  for 
discretizations  derived  from  16-point  and  36-point  prod- 
uct rules,  respectively.  The  accuracy  of  the  results  can  be 
estimated  from  the  fact  that  results  computed  from  the 
different  discretizations  agree  to  better  than  2 digits  in 
the  TM  case  and  better  than  3 digits  in  the  TE  case.  As 
for  the  previous  problem,  however,  the  accuracy  of  the 
propagation  constants  for  the  other  modes  declines  with 
increasing  mode  number. 

Plots  of  the  vector  modes  are  given  in  Figures  1 through 
5.  Small  arrows  indicate  the  local  direction  of  the  (trans- 
verse component  of  the)  electric  field  in  the  aperture.  The 
background  shading  represents  the  corresponding  scalar 
potential. 

Figure  1 shows  the  two  TEM  modes  of  the  guide.  The 
first  was  derived  from  the  solution  to  the  2d  Laplace  equa- 
tion assuming  a unit  potential  on  the  left  interior  con- 
ductor and  zero  potential  on  all  other  boundaries.  The 
second  plot  is  essentially  its  mirror  image.  The  lowest 
eight  TM  modes  are  shown  in  Figures  2 and  3 and  the 
lowest  eight  TE  modes  are  shown  in  Figures  4 and  5. 

2The  number  of  computed  modes  in  this  case  equals  the  number 
of  quadrature  points. 


Fig.  1.  TEM  modes  of  a 5A  x 3A  waveguide  containing  two  interior 
conductors. 

VI.  Conclusions 

Our  method  is  based  on  the  following  observations:  the 
scalar  traveling  modes  of  a waveguide  can  be  obtained 
by  diagonalizing  an  integral  operator  whose  kernel  is  the 
Green  function  for  the  2 d scalar  Helmholtz  equation  in- 
side the  waveguide  aperture;  the  2d  electrostatic  modes 
are  solutions  to  the  2d  Laplace  equation  inside  the  aper- 
ture; and,  the  electromagnetic  modes  can  be  obtained  by 
taking  gradients  of  these  scalar  modes. 

We  showed  how  to  construct  discretized  representa- 
tions of  the  various  integral  operators  using  the  locally 
corrected  Nystrom  method  and  presented  results  from  a 
software  implementation  of  this  method. 

Improvements  are  possible  in  a number  of  areas.  For 
example,  if  we  were  to  set  up  the  2d  Laplace  problem 
using  a double-layer  potential  instead  of  a single-layer 
potential,  the  resulting  integral  equation  would  be  sec- 
ond kind  rather  than  first  kind.  Second  kind  integral 
equations  are  better  conditioned  and  generally  lead  to 
more  accurate  solutions,  especially  for  high  spatial  fre- 
quency modes.  If  we  went  one  step  further  and  combined 
the  double  layer  potential  with  a single  layer  potential 
(i.e.,  employ  a combined  source  formulation),  the  result- 
ing equation  would  be  second  kind  and  would  also  be  in- 
sensitive to  spurious  resonances.  Similar  considerations 
apply  to  the  integral  equations  representing  the  Green 


90 


ACES  JOURNAL,  VOL.  17,  NO.  1,  MARCH  2002,  SI:  APPROACHES  TO  BETTER  ACCURACY/RESOLUTION  IN  CEM 


functions  for  the  scalar  traveling  modes. 


Our  original  objective  was  to  compute  an  accurate,  dis- 
cretized representation  of  the  waveguide  integral  equa- 
tion3 in  order  to  model  waveguide  apertures  and  excita- 
tions in  general  antenna  and  scattering  problems.  This 
requires  accurate  representations  of  the  electromagnetic 
modes  on  the  aperture.  The  results  as  presented  here  are 
not  sufficiently  accurate  for  this  purpose.  Consequently, 
this  paper  should  be  regarded  as  a report  on  a work  in 
progress. 


We  are  grateful  for  support  from  the  Raytheon  Com- 
pany. 
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3 The  waveguide  integral  equation  [3]  is  an  exact  specification  of 
the  relationship  between  the  transverse  electric  and  magnetic  fields 
in  any  cross  section  of  a waveguide.  When  a waveguide  takes  part 
in  a scattering  problem,  the  waveguide  integral  equation  plays  the 
role  of  a (non-local)  boundary  condition  relating  Ej_  and  Hj_  (or 
M =n  x E and  J — n X H). 


Fig.  4.  First  four  TE  modes  of  a 5A  x 3A  waveguide  containing  Fig.  5.  Second  four  TE  modes  of  a 5A  x 3A  waveguide  containing 
two  interior  conductors,  two  interior  conductors. 


